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Abstract 

In a fully 3-D system such as a stellarator, the toroidal mode number n ceases to be a good quantum 
number — all ns within a given mode family being coupled. It is found that the discrete spectrum of unstable ideal 
MHD (magnetohydrodynamic) instabilities ceases to exist unless MHD is modified (regularized) by introducing 
a short-perpendicular-wavelength cutoff. Attempts to use ray tracing to estimate the regularized MHD spectrum 
fail due to the occurrence of chaotic ray trajectories. In quantum chaos theory, strong chaos in the semiclassical 
limit leads to eigenvalue statistics the same as those of a suitable ensemble of random matrices. For instance, the 
probability distribution function for the separation between neighboring eigenvalues is as derived from random 
matrix theory and goes to zero at zero separation. This contrasts with the Poissonian distribution found in 
separable systems, showing that a signature of quantum chaos is level repulsion. In order to determine whether 
eigenvalues of the regularized MHD problem obey the same statistics as those of the Schrodinger equation in 
both the separable 1-D case and the chaotic 3-D cases, we have assembled data sets of ideal MHD eigenvalues 
for a Suydam-unstable cylindrical (1-D) equilibrium using Mathematica and a Mercier- unstable (3-D) equilibrium 
using the CAS3D code. In the 1-D case, we find that the unregularized Suydam-approximation spectrum has an 
anomalous peak at zero eigenvalue separation. On the other hand, regularization by restricting the domain of kx 
recovers the expected Poissonian distribution. In the 3-D case we find strong evidence of level repulsion within 
mode families, but mixing mode families produces Poissonian statistics. 
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1. Introduction 

In ideal MHD the spectrum of the growth rates, 
7, of instabilities is difficult to characterize math- 
ematically because the linearized force operator 
is not compact pQ. This gives rise to the possi- 
bility of a dense set of accumulation points (de- 
scriptively called the "accumulation continuum" 
by Spies and Tataronis 2 though more correctly 
termed [3] the essential spectrum). 

The continuous spectrum in quantum mechan- 



ics arises from the unboundedness of configura- 
tion space, whereas the MHD essential spectrum 
arises from the unboundedness of Fourier space — 
there is no minimum wavelength in ideal MHD. 
This is an unphysical artifact of the ideal MHD 
model because, in reality, low-frequency instabil- 
ities with |kj_| much greater than the inverse of 
the ion Larmor radius, a;, cannot exist (where 
is the projection of the local wavevector into the 
plane perpendicular to the magnetic field B). 



Perhaps the greatest virtue of ideal MHD in fu- 
sion plasma physics is its mathematical tractabil- 
ity as a first-cut model for assessing the stabil- 
ity of proposed fusion-relevant experiments with 
complicated geometries. For this purpose a sub- 
stantial investment in effort has been expended 
on developing numerical matrix eigenvalue pro- 
grams, such as the three-dimensional (3-D) TERP- 
SICHORE |U and CAS3D codes. These solve 
the MHD wave equations for perturbations about 
static equilibria, so that the eigenvalue uo 2 = —j 2 
is real due to the Hermiticity (self-adjointness jH) 
of the linearized force and kinetic energy opera- 
tors. They use finite-element or finite-difference 
methods to convert the infinite-dimensional Hilbert- 
space eigenvalue problem to an approximating finite- 
dimensional matrix problem. 

In order properly to verify the convergence of 
these codes in 3-D geometry it is essential to un- 
derstand the nature of the spectrum — if it is quantum- 
chaotic then convergence of individual eigenvalues 
cannot be expected and a statistical description 
must be used. 

It is the thesis of this paper that the language 
of quantum chaos j7] theory indeed provides such 
a statistical framework for characterizing MHD 
spectra in that it seeks to classify spectra statisti- 
cally by determining whether, and to what degree, 
they belong to various universality classes. 

In the cylindrical case the eigenvalue problem 
is separable into three one-dimensional (1-D) eigen- 
value problems, with radial, poloidal, and toroidal 
(axial) quantum numbers Z, m, and n, respec- 
tively. It is thus to be expected a priori that 
the spectrum will fall within the standard quan- 
tum chaos theory universality class for integrable, 
non-chaotic systems 0. In particular, it is to be 
expected that the probability distribution func- 
tion for the separation of neighboring eigenvalues 
is a Poisson distribution. However, the nature of 
the MHD spectrum is quite different from that 
of the typical quantum, microwave and acoustic 
systems normally dealt with in quantum chaos 
theory and it is necessary to test this conjecture 
by explicit calculation. In fact we find that the 
result depends on the method of regularization. 

We first present the eigenvalue equation for a 
reduced MHD model of a large-aspect-ratio (effec- 
tively cylindrical) stellarator. We study a plasma 
in which the Suydam criterion |H| for the stabil- 



ity of interchange modes is violated, so the num- 
ber of unstable modes tends to infinity as the 
small- wavelength cutoff tends to zero. To com- 
pute large-m eigenvalues we transform to a Schrodinger- 
like form of the radial eigenvalue equation 
which has essentially the same form in config- 
uration (r) space as in Fourier (k r ) space, thus 
allowing easy regularization by restricting the k r 
domain. To simplify even further we approximate 
the effective potential by a parabola, thus yielding 
the quantum harmonic oscillator equation, solv- 
able in parabolic cylinder functions 

Real, finite-aspect-ratio stellarators are fully 3- 
D and their ideal-MHD spectra may be expected 
a priori to fall within the universality class ap- 
propriate to time-reversible quantum chaotic sys- 
tems, where the spectral statistics are found to be 
the same as for a Gaussian orthogonal ensemble of 
random matrices |7| in regions where ray tracing 
reveals chaotic dynamics At the end of this 
paper we give a brief report of 3-D calculations 
peformed with the CAS3D code on a Mercier- 
unstable, high-mirror-ratio, high-iota equilibrium 
representing a Wendelstein 7-X (W7-X) stellara- 
tor variant )12j . 

2. One-dimensional model eigenvalue equa- 
tion 

In this paper we study an effectively circular- 
cylindrical MHD equilibrium, using cylindrical co- 
ordinates such that the magnetic axis coincides 
with the z-axis, made topologically toroidal by 
periodic boundary conditions. Thus z and the 
toroidal angle £ are related through Q = z/R , 
where i?o is the major radius of the toroidal plasma 
being modeled by this cylinder. The poloidal 
angle is the usual geometric cylindrical angle 
and the distance r from the magnetic axis labels 
the magnetic surfaces (the equilibrium field being 
trivially integrable in this case). The plasma edge 
is at r = a. 

In the cylinder there are two ignorable coordi- 
nates, 9 and so the components of £ are com- 
pletely factorizable into products of functions of 
the independent variables separately. In particu- 
lar, we write the r-component as 

r£ r = exp(im6>) exp(— inQip{r) , (1) 

where the periodic boundary conditions quantize 
m and n to integers and we choose to work with 
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the stream function <p(r) = r£ r (r). 

Since the primary motivation of this paper is 
stcllarator physics, we use the reduced MHD or- 
dering for large- aspect stellarators p3E], aver- 
aging over helical ripple to reduce to an equiva- 
lent cylindrical problem |15lllfij . The universality 
class should be insensitive to the precise choice of 
model as long as it exhibits the behavior typi- 
cal of MHD instabilities in a cylindrical plasma, 
specifically the existence of interchange instabili- 
ties and the occurrence of accumulation points at 
finite growth rates. 

Defining A = lo 2 we seek the spectrum of A- 
values satisfying the scalar equation 



Lip = XMip 



(2) 



under the boundary conditions ip(0) = at the 
magnetic axis and tp(l) = 0, appropriate to a per- 
fectly conducting wall at the plasma edge (using 
units such that r = 1 there). 

The operator M = — and L is given by 



L = 
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where G is a Suydam stability parameter (> 1/4 
for instability 8 ), proportional to the pressure 
gradient p'(r) and the average field line curvature 
d 

In this paper we use the notation / = rf'(r) 
for an arbitrary function /, so i = rdb/dr is a 
measure of the magnetic shear and * measures the 
variation of the shear with radius. 

We observe some differences between Eq. (jSJ 
and the standard quantum mechanical eigenvalue 
problem Hip = Eip. One is of course the physi- 
cal interpretation of the eigenvalue — in quantum 
mechanics the eigenvalue E = Tioj is linear in 
the frequency because the Schrodinger equation 
is first order in time, whereas our eigenvalue A is 
quadratic in the frequency because it derives from 
a classical equation of motion. 

Another difference is that Eq. J5J is a gener- 
alized eigenvalue equation because M is not the 
identity operator. This is one reason why it is 
necessary to treat the MHD spectrum explicitly 
rather than simply assume it is in the same uni- 
versality class as standard quantum mechanical 
systems. 
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Figure 1: The rotational transform *(r) = l/q(r) 
with *o = 0.45, *2 = 0.2. All distinct rational 
magnetic surfaces fi = n/m are shown for m up 
to 10. 

Equation J2J is very similar to the normal mode 
equation analyzed in the early work on the inter- 
change growth rate in stellarators by Kulsrud |15| . 
However, unlike this and most other MHD stud- 
ies we are concerned not with finding the highest 
growth rate, but in characterizing the complete 
set of unstable eigenvalues. 

Suydam instabilities occur only for values of 
m and n such that n — mt vanishes. For the 1-D 
numerical work in this paper we use a parabolic 
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transform profile t 
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Given a rational fraction fj. 



*2?" 2 as illustrated in 



n^/m^ in the 



interval [t(0),i(a)] (where n^and m M are mutu- 
ally prime) there is a unique radius r^ 



*(?>) 



such that 
Any pair of integers (m, n) M „ = 



(um^, vn^), v = 1, 2, 3, . . . satisfies the resonance 
condition 



m, 



= 



(4) 



We use a broad pressure profile that is suffi- 
ciently flat near the magnetic axis that the Suy- 
dam instability parameter G goes to zero at the 
magnetic axis, and for which p' vanishes at the 
plasma edge. The resulting G-profile is shown in 

Fig.m 

Defining a scaled radial variable x = m(r — 
r^/r^^we can find the large-m spectrum of Eq. J5J) 
by expanding all quantities in inverse powers of 
to, and equating the LHS to zero order by order. 

In this paper we work only to lowest order in 
1/to, the Suydam approximation. As found by 
Kulsrud |15| . we have the generalized eigenvalue 
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Figure 2: The Suydam criterion parameter G(r) 
(solid line), and the instability threshold 1/4 
(dashed line), showing nearly all the plasma is 
Suydam unstable. 

equation 

£(oyo) = l| f L<°) - W') ^ (0) = , (5) 

TO 2 V / 

where, more explicitly, 
Z"(°) d d 

V = -4(- 2 + r 2 )^ + x 2 + r 2 -G, (6) 

a ax ax 

with r 2 = — X^/i 2 and i and G evaluated at 
Tfj,. Under the boundary conditions <p^> — > as 
7- — > ±oo, Eq. (JSJ can be solved to give a square- 
integrable eigenfunction, with growth rate 7 = tT, 
provided A' ' < is one of the eigenvalues X^j. 
The radial mode number / = 0, 1, 2, . . . denotes 
the number of nodes of the eigenfunction <^( ) = 
¥V,z(0- Note that A^ depends only on ^ = n/m 
and is otherwise independent of the magnitude of 
m and n. 

Restricting attention to unstable modes, so that 
7 = (—A) 1 / 2 is real, we transform Eq. J^J to the 
Schrodinger form [H] 

+ Qfoty = , (7) 

where 

Q = G-\- isech 2 t] - T 2 cosh 2 r\ , (8) 

with r\ defined through x = 7 sinhry/^r^), and 
?/> = (cosh77) 1 / 2 (/3(x). 

From, e.g., Eq. (4.7) of we see that, pro- 
vided the Suydam criterion G > 1/4 is satisfied, 
there is an infinity of 7 eigenvalues accumulating 
exponentially toward the origin from above (so 



the A-values accumulate from below) in the limit 
I —> 00. 

Perhaps less widely appreciated (because m 
and n are normally taken to be fixed) is the fact 
that there is also a point of accumulation of the 
eigenvalues of Eq. (J2J) at each A Mi / as m — > 00 
with I fixed. (Although A^ ^ is infinitely degener- 
ate, we can break this degeneracy by proceeding 
further with the expansion in 1/m, thus showing 
that A M ,2 is an accumulation point.) Since the ra- 
tionals fi are dense on the real line, there is an 
"accumulation continuum" between 7 = and 
the maximum growth rate, 7 = 7 max - 

3.Regularization 

The accumulation points of the ideal MHD spec- 
trum found above are mathematically interesting 
but exist only as a singular limit of equations 
containing more physics, including finite-Larmor- 
radius (FLR) effects and dissipation, that regu- 
larize the spectrum. 

In order to proceed further we need to be ex- 
plicit about the nature of this singular limit. As 
we are primarily concerned with the universality 
class question, we seek only a minimal modifica- 
tion of Eq. J5J that has some physical basis but 
makes as little change to ideal MHD as possible. 
To preserve the Hermitian nature of ideal MHD 
we cannot use the drift correction used for esti- 
mating FLR stabilization of interchange modes 
by Kulsrud ^3- However it is possible to ef- 
fect a pseudo-FLR regularization of ideal MHD 
by restricting to a disk of radius less than the 
inverse ion Larmor radius. In our nondimension- 
alized, large-aspect ratio model this implies 

(k 2 e + k 2 r yi 2 P * < 1 , (9) 

where k r and kg are the radial and poloidal com- 
ponents of the wavevector, respectively, and p* 
is the ion Larmor radius (at a typical energy) in 
units of the minor radius. 

To apply Eq. precisely we need to relate k r 
and kg to the eigenvalue problem discussed above. 
From Eq. we see that kg — m/r. We define 
k r as the Fourier variable conjugate to r. Fourier 
transformation of Eq. (j2J is only practical in the 
large-m limit, when modes are localized near the 
resonant surfaces r = r^, which is why we have 
restricted the discussion to leading order in the 
1/m expansion. 
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Using the stretched radial coordinate x = m(r— 
t^/th we define k r = run/r^, where K is the 
Fourier-space independent variable conjugate to 
x. With the substitutions d/dx i— * i— > id/dn, 
and using the fact that nd/dn and (d/dn)n = 
1 + nd/dn commute, Eq. Q transforms to 



dn ^ 



+ k z )—+T z {\ + k z )-G 
dn 



<p K =0 . 



(10) 

The transformation k = sinh r\ then leads back to 
Eq. (JZJ), with 77 now to be interpreted as a dis- 
torted Fourier-space independent variable, rather 
than as a real-space coordinate! 

Equation implies that Eq. (|10|) is to be 
solved on the domain — K max < k < K max where 



1 



1/2 



This exists provided \m\ < TO max , where 



(11) 



(12) 



Analogously to quantum mechanical box-quantization 
we use Dirichlet boundary conditions at ±K max . 

4. Spectral statistics in the 1-D case 

As only the qualitative nature of the spectrum 
is important, we approximate the function Q by 
Q(0) + iQ"(0)?7 2 , so Eq. can be solved in 
parabolic cylinder functions JOj ■ We find the dis- 
persion relation 



G-r 2 



(4r 2 - i)V2 



(13) 



where v = I in the unregularized case, ft max = 
??max = 00. In the even-^, regularized case v may 
be found by solving for a zero of M 1/2, (4r 2 — 
l) 1 / 2 77 2 nax /2), where M is Kummer's function. For 
I = 0, v becomes exponentially small as 7y max — > 
00, which allows an approximate regularization 
formula to be derived. 

We study the spectrum between the maximum 
I = 1 growth rate, 7 max (Z = 1) = max^^i, 
and the maximum overall growth rate, 7 max = 
max fI 7 (lj o- Only the I = modes exist in this 
range of 7, which corresponds to the range in \x 
between /i m i n « 0.522 and /i max ~ 0.628. Through- 
out this range T is > 1/2, so that Q has a single 
minimum !) and the quadratic approximation of 




Figure 3: The histogram shows an estimate, based 
on a data set of about 32, 000 unregularized eigen- 
values, of the probability distribution function for 
the eigenvalue separation s. The plot is domi- 
nated by the spike at s = 0. 

this section is appropriate. In this range there are 
only four low-order rationals n/m with m < 10. 

Taking p» = 0.001, all pairs of integer values to, 
n in the fan-shaped region 1 < to < TO max (n/TO.), 
Mmin < n/m < /i max were evaluated, giving an 
initial dataset of over 32,000 points (to, n). The 
corresponding set of unregularized eigenvalues was 
calculated by solving Eq. (|13JI with v = and the 
eigenvalues were sorted and numbered from the 
top to give the integrated density of states "stair- 
case" function N(j). 

The curve 0.3523-9.5733 x 10- n iV 2 - 1.1625 x 
10 -20 Af 4 was found to give a good fit to the smoo- 
thed behavior of this function. Inverting this func- 
tion gives the smoothed function ^(7) which is 
used to "unfold" spectra by defining a new 
"energy eigenvalue" E = N(j), such that N(E) 
increases linearly on average. 

This means that the average separation of eigen- 
values is now unity, making comparison with spec- 
tra from other physical systems meaningful and 
allowing universal behavior to become apparent if 
present. However, Fig.[3]shows that the probabil- 
ity distribution of eigenvalue spacings s is far from 
universal for the unregularized Suydam spectrum, 
exhibiting a dclta-function-like spike at s = 0. 
This is presumably because, although we have 
truncated the spectrum in to, we have not re- 
moved the degeneracies arising for low-order ra- 
tionals /j in the range /i m i n < /1 < /i m ax- 

Figure 0] on the other hand shows that when 
a similar procedure is applied to the regularized 
spectrum (retaining only regularized eigenvalues 
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Figure 4: Eigenvalue spacing distribution for a 
data set of about 12, 000 regularized eigenvalues. 
The exponential curve shows the Poisson distri- 
bution. 




Poisson process 
chaotic process 
N=0 family 



Poisson process 

chaotic process 

N=2 family 




Poisson process 
chaotic process 
N=1 family 




Poisson process 

chaotic process 

ALL families 



above 7 max (^ = 1), the universal Poisson distri- 
bution expected from a separable system is ob- 
tained to a good approximation, thus leading to 
the expectation that generic quantum chaos the- 
ory is applicable once any physically reasonable 
regularization is performed. 

Further support for this hypothesis is obtained 
from a CAS3D study of a W7-X variant equilib- 
rium with a nonmonotonic, low-shear transform 
profile (* axis = 1.1066, t min = 1.0491, i odg c = 
1.0754). As seen from Fig. O when the statistics 
are analyzed within the three mode families the 
eigenvalue spacing distribution function is closer 
to the Wigner conjecture form found for generic 
chaotic systems [J] than to the Poisson distribu- 
tion for separable systems, as might be expected 
from 11 . However, when the spectra from the 
three uncoupled mode familes are combined, there 
are enough accidental degeneracies that the spac- 
ing distribution becomes close to Poissonian. 
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